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ABSTRACT 

For the next generation of radio interferometric telescopes it is of paramount importance to 
incorporate wide field-of-view (WFOV) considerations in interferometric imaging, otherwise 
the fidelity of reconstructed images will suffer greatly. We extend compressed sensing tech- 
niques for interferometric imaging to a WFOV and recover images in the spherical coordinate 
space in which they naturally live, eliminating any distorting projection. The effectiveness 
of the spread spectrum phenomenon, highlighted recently by one of the authors, is enhanced 
when going to a WFOV, while sparsity is promoted by recovering images directly on the 
sphere. Both of these properties act to improve the quality of reconstructed interferomet- 
ric images. We quantify the performance of compressed sensing reconstruction techniques 
through simulations, highlighting the superior reconstruction quality achieved by recovering 
interferometric images directly on the sphere rather than the plane. 
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1 INTRODUCTION 

Incorporating wide field-of-view (WFOV) contributions in the im- 
ages reconstructed from radio interferometric observations is be- 
coming increasingly important. Next-generation radio interferom- 
eters, such as the Square Kilometre Arra)[](SKA) ( Carilli & Rawl- 
|ings|2004] >, will inherently observe very large fields of view about 
the pointing direction of the telescope. Wide fields introduce two 
important distinctions to standard interferometric imaging: firstly, 
interferometric images are inherently spherical and planar pro- 
jections necessarily introduce distortions; and, secondly, non-zero 
baseline components in the pointing direction of the telescope must 
be taken into account. If these contributions are ignored, the fidelity 
of reconstructed images will suffer greatly. 

WFOV contributions have been considered by |McEwen &| 
Scaife (2008) in simulating the visibilities observed by an inter- 
ferometer. Full- sky interferometric formalisms were derived using 
a number of different signal representations, including represen- 
tations in real, spherical harmonic and wavelet spaces. Real and 
spherical harmonic space representations were shown to be nu- 
merically infeasible for simulating realistic observations, while a 
fast wavelet space method was developed, reducing the compu- 
tational cost considerably and rendering realistic simulations fea- 
sible. However, the forward problem only was considered in this 
work. Very recently |Carozzi & Woan| ( |2009| > developed a gen- 
eralised radio interferometric measurement equation that is valid 
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for partially polarised sources over a WFOV, however this frame- 
work has not yet been applied in practice. The inverse WFOV 
imaging problem has traditionally been tackled by faceting the sky 
into a number of regions which are sufficiently small that standard 
Fourier imaging is possible (Cornwell & Perley 1992). More re- 
cently, the w-projection algorithm has been developed by Cornwell 
et al. (2008) to incorporate WFOV contributions by applying the 
modulating so-called w-term that appears in the interferometric in- 
tegral as a convolution in the Fourier plane. This provides an order 
of magnitude speed enhancement over traditional facet-based ap- 
proaches. Nevertheless, both of these methods recover planar im- 
ages in the space of directional cosines; this necessarily distorts the 
image. None of the current methods recover images in the spheri- 
cal coordinate space in which they live. In this article we recover 
spherical interferometric images parameterised by colatitude and 
longitude on the celestial sphere. By recovering images defined di- 
rectly on the sphere, and thereby eliminating the distortion due to 
a projection to the plane, the performance of reconstruction is en- 
hanced. We develop reconstruction algorithms in the context of the 
theory of compressed sensing. 

Compressed sensing {Candes et~ l. 2006a b, Candes 2006, 
Donoho 2006, Baraniuk 2007} is a recent development in the 
field of information theory, which goes beyond the usual Nyquist- 
Shannon sampling theorem. It relies on the fact than many sig- 
nals in Nature are sparse (or approximately so) and may be repre- 
sented in a basis requiring many fewer non-zero coefficients than 
the dimensionality of the signal itself. Compressed sensing the- 
ory shows that a sparse signal may be recovered from many fewer 
measurements than Nyquist-Shannon sampling would suggest and 
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thus aims to merge data acquisition and compression. These results 
also hold for signals that are approximately sparse only (i.e. sig- 
nals which contain many coefficients of small but non-zero value), 
so-called compressible signals. 

The first application of compressed sensing to radio interfer- 
ometry was performed by Wiaux et al. ( 2009a ), where the problem 
of image reconstruction from incomplete visibility measurements 
was considered. |Wiaux et al.| < |2009a| > demonstrated the versatility 
of the approach, through the ability to incorporate additional sig- 
nal priors easily, and its superiority relative to standard interfero- 
metric imaging techniques, such as CLEAN (Hogbom 1974). The 
spread spectrum phenomenon, which arises by partially relaxing 
the small field-of-view (FOV) assumption and including a first or- 
der w-term, was introduced by Wiaux et al. (2009b), enhancing the 
performance of image reconstruction for sparsity bases that are not 
maximally incoherent with the measurement basis (see Sec. |2.2| for 
a review of sparsity and measurement bases and their coherence). 
Furthermore, a compressed sensing based approach was developed 
and evaluated by | Wiaux et al.| ( [20T0| to recover the signal induced 
by cosmic strings in the cosmic microwave background, exploit- 
ing the sparse nature of line-like discontinuities due to the string 
signal. All of these works consider uniformly random and discrete 
visibility coverage in order to remain as close to the theory of com- 
pressed sensing as possible. First steps towards more realistic visi- 
bility coverages have been taken by Suksmono (2009) and Wenger 
|et al.| ( |2010| ), who consider coverages due to specific interferom- 
eter configurations but which remain discrete. These preliminary 
works suggest that the performance of compressed sensing recon- 
structions is unlikely to deteriorate significantly for more realistic 
visibility coverages. 

In this article we generalise the compressed sensing imaging 
techniques developed by Wia ux et al.| ( |2009a| ) and | Wiaux et al.| 
( 2009b ) to a WFOV. To do so we recover interferometric images de- 
fined directly on the sphere, rather than a tangent plane. Recovering 
images in the space in which the signal naturally lives eliminates 
any distortion due to projection. Furthermore, projection effects 
also act to hamper the sparsity of the signal (i.e. act to make it less 
sparse), impeding the performance of compressed sensing based re- 
construction. To remain close to the theoretical compressed sensing 
framework we follow | Wiaux et al.| ( |2009a} and |Wiaux et al.| (|2009b ) 
by considering uniformly random and discrete planar visibility cov- 
erage. We also assume non-zero but constant baseline components 
in the pointing direction of the telescope (i.e. non-zero but constant 
w, where w is defined explicitly in Sec. |2.1| ), in order to study the 
enhancement of reconstruction quality due to the spread spectrum 
phenomenon. This assumption allows us to discard considerations 
related to specific interferometer configurations and to study the 
impact of the spread spectrum phenomenon at light computational 
load. As we shall see, the effectiveness of the spread spectrum phe- 
nomenon is improved in the WFOV setting considered. Due to all 
of these considerations, the quality of reconstruction on the sphere 
is enhanced considerably relative to planar reconstructions. 

The remainder of this article is structured as follows. In Sec. [2] 
we review radio interferometric imaging in the context of com- 
pressed sensing. In Sec. [5] we generalise these techniques to a 
WFOV. Simulations are presented in Sec. [4] to evaluate the perfor- 
mance of the compressed sensing based reconstruction of spherical 
images compared to planar images. Concluding remarks are made 
in Sec. [5] 



2 PLANAR INTERFEROMETRIC IMAGING 

In this section we review radio interferometric imaging in the con- 
text of compressed sensing. Radio interferometry is reviewed in 
the WFOV setting, before small FOV assumptions are discussed 
explicitly. A brief review of compressed sensing is given, highlight- 
ing various reconstruction techniques and the importance of spar- 
sity and coherence on reconstruction performance. Finally, these 
frameworks are merged and compressed sensing techniques for ra- 
dio interferometric imaging are discussed. 

2.1 Radio interferometry 

In order to image a region on the sky, all radio telescopes of an in- 
terferometric array are oriented towards the same pointing direction 
s on the unit celestial sphere S 2 . The FOV observed by the interfer- 
ometer is limited by the primary beam of the telescope A(r), where 
the beam is defined relative to the pointing direction, i.e. r = s - s 
for arbitrary directions s. The sky intensity to be imaged x(r) is 
defined in the same coordinates. The complex visibility measured 
by each telescope pair of the interferometer is given by ( Thompson 
|et al.|2001) 

y(b)= f A(t)x(t)z-' M t dn(r), (1) 

where dQ(s) = sin 6 d6 dtp is the usual rotation invariant mea- 
sure on the sphere and (6, if) denote the spherical coordinates of 
s, with colatitude 6 e [0, n] and longitude <p e [0, 2tt). The mea- 
sured visibility depends only on the relative positions between tele- 
scope pairs, denoted in units of wavelengths by the baseline vector 
b = (u, v,w). As the Earth rotates relative to the celestial sphere 
the interferometer tracks the source position as it traverses the sky, 
with each position corresponding to a rotation of the baseline. In 
this manner, visibility measurements are accumulated for various 
baselines, with each antenna pair generating an elliptical track of 
baseline values over the course of the observation. Visibility cover- 
age is therefore incomplete. 

The beam and sky intensity are typically expressed in the same 
Cartesian coordinate system as the baseline, which is centred on the 
Earth and aligned with s . Consider the coordinates (/, m, n) of s in 
this system, noting n = n(l) = (1 - ||/||2) 1/2 with / = (l,m) and 
||/|| 2 = I 2 + m 2 . The primary beam and sky intensity may then be 
seen as functions of /. Following this change of coordinates the 
visibility integral of (T} reads 

y(u,w)= f A(l)x(l)e- i2 « [u - l+w ^ l) - 1)] , (2) 
Jd2 n(l) 

where u = (u, w) and we have noted that r = (/, n(l) - 1). Through 
the change of coordinates the invariant measure on the sphere be- 
comes d£l(s) = ft _1 (7)d 2 /, where d 2 / = d/dra is the canonical 
invariant measure on the plane, and the integration is now per- 
formed over the unit disk D 2 . Note that (2) remains general and 
does not rely on any small FOV assumption. However, the direc- 
tional cosines / = (/, m) express a projection onto the equatorial 
plane of the signal, which is inherently defined on the celestial 
sphere. This projection is trivial in the continuous setting but will 
become important once we reach the discrete setting required for 
interferometric imaging. 

Often small FOV assumptions are made, in which case it is 
convenient to represent the so-called w-term 

c(\\i\\ 2 ) = e 2 ^-^^) (3) 
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explicitly, from which it follows 

y(u,w)= f A(l)x(l)C(\\l\\ 2 )^ u l ^ (4) 
Jd2 nil) 

Two types of approximation regarding small FOVs are used to 
relate ([4} to the Fourier transform of the beam-modulated inten- 
sity A(l)x(l). Firstly, the assumption \\l\\^ <§c 1 implies n(l) ^ 1, 
so that the Jacobian term n~ l (l) is reduced to unity. Alternatively, 
n~ l (l) may be absorbed into the beam. Secondly, various assump- 
tions are made to approximate the w-term. A zeroth order approx- 
imation of the w-term is made by assuming w «: 1, so that 
COI/H2) - 1- When incorporating both approximations, ([4} reduces 
to the Fourier transform and standard Fourier imaging may be 
used to recover images from measured visibilities. Alternatively, 
a first order approximation Hfll* w «: 1 reduces the w-term to 
C(||/|| 2 ) - e i7rvv||/|l 2 = Ci(||/|| 2 ). This assumption gives rise to the 
linear chirp modulation responsible for the spread spectrum phe- 
nomenon (Wia ux et al.|2009b| . Although not considered here, di- 
rection dependent beam effects may also introduce a phase mod- 
ulation, which would provide an alternative source of the spread 
spectrum phenomenon. Since we intend to consider WFOVs, no 
small FOV assumptions will be made herein; we include the Jaco- 
bian term n~ l (I) and consider the full w-term given by (p). 



2.2 Compressed sensing 

Compressed sensing ( Candes et al. 2006a b ; Candes 2006 ; Donoho 
2006, Baraniuk 2007 ) is concerned with the recovery of sparse or 
compressible signals from a small number of measurements. The 
signal to be recovered is defined by its Nyquist-Shannon sampling, 
denoted by the vector x e C N , and is assumed to be sparse in some 
orthogonal basis {^ f }i<i<#, where e C N , V/. The signal may then 
be represented by its coefficients a e C N in this basis: x = ^a, 
where *F e £ NxN is the N x N matrix with columns Formally, x 
is said to be sparse or compressible if a contains K <z: N non-zero 
or significant elements respectively. Hereafter, we refer to enhanc- 
ing (hampering) sparsity as synonymous with decreasing (increas- 
ing) the sparsity value K. Measurement of x is assumed to be made 
by the projection onto measurement basis vectors {0 r h< r <M> 

for 

M measurements, belonging to an orthogonal sensing basis where 
<j) r eC N , Vr. This is a very flexible framework which allows a wide 
range of acquisition procedures to be modelled. The vector of mea- 
surements y e C M may then be expressed by 

y = Ox + n = O^or + n , (5) 

where O e £ MxN is the M x N matrix withs rows <p r and n e C M 
represents noise. Compressed sensing suggests that x can be re- 
covered with a number of measurements M ~ K «: N. However, 
recovering x in this setting involves solving the inverse problem 
which becomes ill-posed for M < N. 

Compressed sensing techniques are generally based on global 
minimisation problems, which are solved by greedy algorithms or 
convex non-linear optimisation algorithms. The ill-posed inverse 
problem described above can be defined by a constrained optimisa- 
tion problem explicitly regularised by a sparsity or compressibility 
prior. This results in the Basis Pursuit denoising (BP) problerrj^] 
which consists of minimising the i\ -norm of the coefficients of the 
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signal in the sparsity basis ||or||i under a constraint on the ^ 2 -norm 
of the residual noise \\y - O^or^: 

min |M|i such that \\y - < e . (6) 

a 

Recall that the ^-norm is simply given by the sum of the abso- 
lute values of the elements of a vector = Yli=\ whereas 
the £2 -norm is the standard norm defined previously. The constraint 
6 may be related to a residual noise level estimator. Assuming in- 
dependent identically distributed Gaussian noise, a residual noise 
level estimator is given by twice the negative logarithm of the like- 
lihood associated with the candidate reconstruction, which follows 
a ^-distribution with 2M degrees of freedom. The measurement 
constraint € may then be chosen to correspo nd to some (100p)th 
percentile of the ^-distribution, for < p < 1 (Wi aux et al.|20 09a). 
The BP problem may be solved by the application of non-linear, it- 
erative convex optimisation algorithms (e.g.Combettes & Pesquet 
|2007| ). If the solution of the optimisation problem is denoted or*, 
then the signal is reconstructed through the synthesis x* BP = 

The BP denoising problem is appropriate for signals sparse 
or compressible in a given basis. However, many signals in Nature 
are also sparse or compressible in the magnitude of their gradient, 
in which case the Total Variation (TV) minimisation problem ap- 
plies ( Candes et al. 2006a ). The TV problem involves replacing the 
^i-norm sparsity prior in the BP problem with a prior on the TV 
norm of the signal itself: 

min IMItv such that \\y - <bx\\ 2 < e , (7) 

where the TV norm is defined by the ^-norm of the gradient of the 
signal ||x||tv = ||Vjc||i ( |Rudin et al.|1992l|Chambolle|2004| ). The TV 
problem may also be solved by the application of non-linear, itera- 
tive convex optimisation algorithms (e.g. Chambolle 2004, Durand 
|et al.|2010| >. The signal is directly recovered from the solution to 
the optimisation problem, denoted x* TW . 

Finally, we review the two fundamental criteria that drive the 
performance of compressed sensing reconstructions. The first has 
already been central to our discussion of compressed sensing: spar- 
sity. The more sparse a signal the fewer measurements required to 
recover it, or similarly, the better the reconstruction quality for a 
given number of measurements. In addition, the matrix = 
must satisfy the restricted isometry property (RIP) ( Candes et al.| 
2006a b, Candes 2006 ) to ensure accurate recovery. It has been 
shown that this property can be satisfied by acquiring random mea- 
surements (i.e. random visibility coverage in the terminology of 
interferometry), if the measurement and sparsity bases are inco- 
herent. Coherence is therefore the second criterion driving recon- 
struction performance; as the coherence between the two bases in- 
creases, the reconstruction performance degrades. Incoherence en- 
sures that the rows of O cannot sparsely represent the columns of 
X F, ensuring that signal content is sufficiently probed by random 
measurements. Note that coherence is only defined strictly in the 
presence of a sparsity basis, hence it cannot be studied for the TV 
problem (nevertheless, approximate coherence arguments may still 
be made in this setting to gain intuition). The mutual coherence be- 
tween the measurement and sparsity bases is defined by the maxi- 
mum absolute inner product of all combinations of their normalised 
basis vectors: 

MO,¥) = max|<0 r |^.>|, 

with the incoherence given by the reciprocal /z -1 (0,*F). Note that 
the Dirac and Fourier bases are maximally incoherent. The number 
of measurements required to satisfy the RIP and ensure accurate 
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recovery, satisfies the bound 

M > c K // 2 (0,¥) AHog 4 N , 

for some constant c (note that the bound is not tight). Although 
useful for theoretical and intuitive considerations, the mutual co- 
herence is a blunt numerical instrument as it captures only the most 
extreme inner product between the measurement and sparsity bases 
(Tropp 2004 ). An alternative measure of coherence is given by the 
cumulative coherence (Romberg |2009| > 



v(<D,«P,n = mffl(£|<*rl*i>| 

T \ ieT 



1/2 



where T represents the set of all sparsity basis vectors that con- 
tribute to the signal of interest. The cumulative coherence is there- 
fore signal dependent and incorporates sparsity information, which 
renders it more robust than mutual coherence to departures from 
the pure compressed sensing framework. The number of measure- 
ments required for signal recovery may also be expressed in terms 
of cumulative coherence and can be shown to evolve as ( Romberg 
|2009l|Rauhut|20T0t 



(8) 



2.3 Interferometric imaging 

We express the interferometric framework discussed in Sec. |2.1| 
in a discrete setting which is amenable to the compressed sens- 
ing reconstruction techniques discussed in Sec. |2.2[ follow ing the 
approach taken by |Wiaux eFaL] < |2009a| ) and |Wiaux et al.| ([2009b ) 
previously. Although we consider WFOV considerations we restrict 
ourselves to baselines with constant w, thus restricting visibilities 
to a Fourier planeQ In reality, w will vary and the impact of the 
spread spectrum phenomenon will lie somewhere between the ex- 
treme cases that we consider here. Nevertheless, this assumption 
allows us to probe the expected impact of the spread spectrum 
phenomenon at light computational load and to discard consid- 
erations related to specific interferometer configurations. Consid- 
ering the signal projected onto the equatorial plane, we consider 
the / = (/, m) coordinate system discretised on a uniform grid of 
N = N l/2 xN 1/2 points l t e R 2 in real space, with integer 1 < i < N, 
and the corresponding grid of discrete spatial frequencies u t defined 
by Nyquist-Shannon sampling. The band-limited intensity and pri- 
mary beam functions are defined on the spatial grid and denoted 
by x, A e C N , respectively. The complex w-term is also defined on 
this grid C e C N . Since the w-term is complex, the Fourier trans- 
form of C(\\l\\ 2 )A(l)x(l) does not bear any specific symmetry in the 
Fourier plane, even for a real beam and real source intensity. Fol- 
lowing W iaux etaT] ( [2009a] ) and | Wiaux et al.| ( [2009b| ), we assume 
that the spatial frequencies u probed by the interferometer fall on 
the discrete grid points u t . The Fourier coverage provided by the 
M spatial frequencies probed u r , with integer 1 < r < M, can be 
identified by a binary mask in the Fourier plane equal to unity for 
each spatial frequency probed and zero otherwiseJ^The visibilities 



3 The w-projection algorithm JCornwell et al.|2008) could be applied in 
future to remove this restriction, although the extension of the spread spec- 
trum phenomenon to a var ying w would need to be studied extensively. 

Conventional gridding (Thompson et al. 12001) of the visibilities mea- 
sured at continuous spatial frequencies would result in a mask with non- 
binary weights, which could be incorporated easily in the framework de- 
scribed here. 



measured are denoted by the vector y e C M , which may be affected 
by complex noise n e C M . 

The visibility integral ^ may be represented in this discrete 
setting by the linear system 



y = O p Xp + n , 



(9) 



with 



WMFCA, 



where Gaussian noise and a whitening operator have also been in- 
cluded. The measurement matrix O p e £ MxN identifies the com- 
plete linear relation between the signal and the visibilities. The ma- 
trix A e £ NxN i s the diagonal matrix implementing the primary 
beam, with the Jacobian term n~ l {l) due to the change to Carte- 
sian coordinates absorbed. The matrix C e £ NxN is the diagonal 
matrix implementing the modulation by the w-term. The unitary 
matrix F e ( 
matrix M e 



implements the discrete Fourier transform. The 
Y is the rectangular binary matrix implementing 
the mask that characterises visibility coverage, containing one non- 
zero value on each row only, at the index of the Fourier coefficient 
corresponding to the spatial frequency probed by each interfero- 
metric measurement. We also augment the measurement matrix by 
a whitening operator W e R MxM ( [Wiaux et al.|2009a| ). Whitening 
corresponds to dividing each measured visibility by the standard 
deviation of the corresponding noise component so that the final 
observed visibilities are then affected by independent identically 
distributed Gaussian noise with unit variance. The subscript p is 
used to denote the planar setting since we subsequently generalise 
this framework to functions defined on the sphere (throughout sub- 
scripts p and s are used to denote the plane and sphere respectively). 
For incomplete visibility coverage M < N, (5} defines an ill-posed 
inverse problem, which we solve using the BP and TV reconstruc- 
tion approaches described in Sec. |2.2| Finally, let us note that the 
compressed sensing framework is defined strictly for orthogonal 
sparsity and sensing bases. The application of the w-term modula- 
tion necessitates an upsampling operation in practice to ensure that 
the modulated signal is Nyquist-Shannon sampled, which breaks 
the orthogonality of the sensing basis. Compressed sensing tech- 
niques for interferometric imaging therefore depart from the theo- 
retical compressed sensing framework but nevertheless have been 
demonstrated to perform very well. 

We conclude this section with a review of the spread spectrum 
phenomenon. As discussed in Sec. |2.2[ the coherence between the 
measurement and sparsity bases is a key criterion effecting the per- 
formance of compressed sensing based reconstructions. In the in- 
terferometric setting described above the measurement basis can 
essentially be identified with the Fourier basis, which will aid our 
intuitive understanding of the spread spectrum phenomenon. In this 
case, the mutual coherence is given by the maximum modulus of 
the Fourier coefficient of the sparsity basis vectors. Consequently, 
an operation that acts to reduce the maximum Fourier coefficient, 
reduces the coherence and thus improves the quality of compressed 
sensing reconstructions. Modulation by the w-term corresponds to 
a norm-preserving convolution in the Fourier plane, spreading the 
spectrum of the sparsity basis vectors, and achieves exactly this. 
Hence, the greater the frequency content of the w-term, the more 
effective the spread spectrum phenomenon. 
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3 SPHERICAL INTERFEROMETRIC IMAGING 

In this section we generalise the compressed sensing imaging tech- 
niques developed by |Wiaux et aT] ( |2009a| ) and |Wiaux et"aT] ( |2009b| ) 
to a WFOV. We recover interferometric images defined directly on 
the sphere, rather than the equatorial plane, which enhances the ef- 
fectiveness of compressed sensing reconstructions. 

The measurement operator transforming the sky intensity de- 
fined on the sphere to visibilities, consists of augmenting the usual 
interferometric measurement operator with an initial projection 
from the sphere to the plane. This initial projection corresponds 
to a change from spherical to Cartesian coordinates in much the 
same way ^ is defined from {T}. Similarly, the framework re- 
mains general and does not rely on any smal FOV assumptions. 
Practically, however, the projection which implements the change 
of variable is complicated by the discrete setting and the desire to 
recover a regular grid on the plane to allow the use of fast Fourier 
transforms (FFTs)[^]ln order to ensure information is not lost, we 
define the resolution of the planar grid so as to support the maxi- 
mum projected frequency content of a band-limited signal on the 
sphere. Although our WFOV framework still involves a projection 
to the plane, this is included in the measurement operator, and so 
will be regularised when solving the interferometric inverse prob- 
lem. Moreover, we recover the sky intensity directly on the sphere, 
where is it most sparse, and it is the signal space that is important 
for the impact of sparsity on the performance of compressed sens- 
ing reconstructions. 

In the remainder of this section we first discuss considerations 
relating to band-limited signals on the sphere and the plane, opera- 
tors used to project the sphere to the plane in a discrete setting and 
discrete gradient operators defined on the sphere (required for TV 
reconstructions on the sphere). Finally, we define interferometric 
imaging in the WFOV setting, while commenting on the impact of 
sparsity, coherence and the spread spectrum phenomenon. 



3.1 Band-limited signals 

We require a regular grid on the plane to enable the application of 
FFTs to reduce considerably the computational load of reconstruct- 
ing images from visibility measurements. To ensure information is 
not lost when projecting a signal from the sphere to the plane in this 
discrete setting, we define a planar grid that supports a band-limit 
corresponding to the maximum projected frequency content of a 
band-limited signal defined on the sphere. The maximum projected 
frequency content arises at the extents of the FOV, where signal 
content defined on the surface of the sphere is projected onto much 
higher frequency content in the / = (/, m) plane. In the typical small 
FOV setting the relationship between the band-limit of a signal on 
the sphere £ max and its tangent plane w max is given by ^ max = 2^w max . 
In the WFOV setting simple geometric considerations at the extent 
of the FOV lead to the relationship 

Cax = 2tt cos(# fov /2) u max , (10) 

where # F ov denotes the angular opening of the FOV, corresponding 
to a planar FOV of L = 2 sin(# FO v/2). The band-limit relation is 
now dependent on the size of the FOV. Note that for a WFOV a 



5 Fast algorithms have been developed to compute a discrete Fourier trans- 
form on non-equispaced spatial frequencies ( Potts et al. 2001), which could 
in principle be used to avoid explicit gridding. 



higher band-limit on the plane is required to support a given band- 
limit on the sphere, as expected intuitively, and as # FO v — > the 
usual relationship for a small FOV is recovered. 

Once the band-limit of a signal is defined, on both the sphere 
and the plane, sampling considerations then dictate the resolution 
of the discrete grid required to accurately represent the signal. Con- 
sequently, once the band-limit of the signal on the sphere and the 
FOV is set, the required band-limit on the plane is determined 
through {TO} , and the sampling resolutions of both the sphere and 
the plane follow. The Nyquist-Shannon sampling theorem on the 
plane states that N p 1/2 = 2w max L, where N p is the number of samples 
on the planar grid. The relationship between the number of samples 
within the FOV on the sphere N s and the harmonic band-limit ^ max , 
depends on the pixelisation of the s phere adopted. We choose the 
HEALPi^pixelisation of the sphere ( |G6rski et aL]2 005) due to the 
equal area of each pixel element and its ubiquitous use in the astro- 
physical community. The resolution of the HEALPix pixelisation is 
controlled by the parameter Af S ide, where the number of pixels at a 
given resolution is specified by Af pix = \2N S i de 2 . Functions that are 
band-limited on the sphere at ^ max can be resolved on a HEALPix 
pixelisation at a resoluti on corresponding to £ max = 3N side , albeit 
with integration error^]( Hivo n et al.|2010| ) (since no exact quadra- 
ture rule exists for HEALPix). For the correspondence ^ max = 2A^ side 
integration errors are minimal ( |Hivon et al.||2 010). Consequently, 
harmonic transforms are typically performed for ^ max = A:A^ side , with 
k ~ 2-3. In order to ensure that we do not favour either the plane or 
sphere in our analysis, we choose k such that in the limit of a small 
FOV, the number of samples on the plane and within the FOV on 
the sphere are the same. It can be shown trivially that 

A/p _ 2k 2 tan 2 Q9 FOV /2) 
N s ~ 3/r 2 (l - cos(# FOV /2)) ' 

where ^ max drops out of the expression. Consequently, we choose 
the value k = V3^/2 to ensure 

lim N p /N s = 1 . 

The ratio N p /N s is plotted in Fig.[T]for various values of k. Notice 
how the ratio increases rapidly with the FOV. Naively, one might 
expect the ratio of sparsities between the plane and sphere K p /K s 
to evolve in a similar manner, which would highlight the superiority 
of spherical reconstructions compared to planar ones (this is neces- 
sarily a very approximate prediction and the exact ratio of sparsities 
will depend highly on the signal examined). 

3.2 Projection operators 

Now that discrete grids are defined on the sphere and the plane, it is 
necessary to determine how to project a signal between these grids. 
In the continuous setting this is trivial and simply corresponds to a 
change of variable, however the matter is complicated in practice 
due to the discrete nature of the sampled signals. 

In order to project onto a regular grid on the plane, it is neces- 
sary to re-grid the pixelisation on the sphere to recover sample val- 
ues at spherical positions that project directly onto the planar grid, 
as illustrated in Fig. [2] We perform a convolution on the sphere 
to achieve this re-gridding. This convolutional re-gridding on the 
sphere is similar to the re-gridding often performed when map- 
ping the visibilities observed at continuous coordinates to a regular 

6 |http : //healpix ■ jpl ■ nasa . gov/| 

7 Integration errors may be reduced by running the analysis iteratively. 
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Figure 1. Ratio of number of samples on the plane to the sphere (N p /N s ) to 
sufficiently sample a band-limited signal defined on the sphere when pro- 
jected onto the plane (#fov is specified in degrees). Curves are plotted for 
a HEALPix pixelisation of the sphere, with £ max = kN s id e , for values k = 3 
(blue/dashed line); k = V3^/2 (black/solid line); k = 2 (red/dot- dashed 
line). To ensure a ratio of unity is obtained as #fov tends to zero, we select 
k = V3 n/2 throughout this work. 



grid, also to afford the use of FFTs (Thomps on et al.|2 001 ). The 
interferometric measurement operator on the sphere is therefore 
augmented by prepending a projection operator P that includes a 
spherical convolution followed by a project from the sphere to the 
tangent plane / = (/, m). 

We consider three convolution kernels on the sphere: a box 
kernel, corresponding to a nearest-neighbour interpolation; a sine- 
like kernel; and a Gaussian kernel. The box kernel is well localised 
on the sphere but has infinite support in harmonic space. The sinc- 
like kernel corresponds to evaluating the spherical harmonic trans- 
form of the sampled signal on the sphere, followed by evaluating 
the signal value at the new sample position from its spherical har- 
monic coefficients (we refer to this as a sinc-like kernel due to anal- 
ogy with the plane). This procedure results in a kernel that is well 
localised in harmonic space but has support in real space that en- 
tends over the entire sphere. These two kernels represent the ex- 
tremes between spatial and harmonic space localisation and each 
produces ringing in the projected signal in the domain in which the 
kernel is not well localised. We therefore seek a kernel that pro- 
vides a compromise between this trade-off and is reasonably well 
localised in both space and frequency. In this article we settle on 
a simple Gaussian kernel, however alternative kernels such as the 
spherical equivalent of the Gaussian-sinc or spheroidal functions 
could also be considered (Thom pson et al.|2001| ). 



n A 




Figure 2. Projection of a sampled signal from the sphere to the plane. In 
order to project onto a regular grid on the plane (to reduce significantly the 
computational load of subsequent analyses through the use of FFTs), a re- 
gridding operation is required. The value of the point on the plane (blue/dark 
dot) is determined by convolving the points on the sphere (red/grey dots) 
with a suitable kernel. 

on the sphere / reads 

which may be computed in harmonic space ( Wandelt et al. 1998 ) 
to eliminate pixelisation concerns. However, such an approach re- 
quires a spherical harmonic transform, which necessarily requires 
global support on the sphere and hence is non-local in nature, and 
also is of high computational cost. 

We define a discrete gradient operator on the sphere by ana- 
logue with the continuous gradient but computed through finite dif- 
ferences. Convolutional re-gridding is performed to obtain samples 
on the sphere at the positions required to compute finite differ- 
ence values. The same Gaussian convolution kernel that is used 
to project the sphere to the plane (as discussed in Sec. |3.2| ) is used. 
Such an approach actually corresponds to computing the gradient 
of a smoothed version of the original signal. However, the smooth- 
ing is minimal in practice and numerical experiments have shown 
that the discrete gradient defined in this manner is a good approxi- 
mation to the continuous gradient computed in harmonic space but 
does not suffer from the global nature and high computational cost 
of the spherical harmonic transform. 



3.3 Gradient operators on the sphere 

A discrete gradient operator must be defined on the sphere in or- 
der to compute the magnitude of the gradient of a function to 
solve the TV minimisation problem of {7} on the sphere. The dis- 
crete gradient operator on the plane is defined simply through finite 
differences ( |Rudin et aL][T99ll |Chambolle|[2004] >. However, it is 
not possible to define a discrete spherical gradient operator on the 
HEALPix pixelisation in this manner since sample positions do not 
lie on both rings of constant latitude and rings of constant longitude 
(only equiangular pixelisations of the sphere satisfy this property). 
One alternative is to consider the continuous gradient operator on 
the sphere, for which the magnitude of the gradient of a function 



3.4 Interferometric imaging 

We are now in a position to define explicitly the WFOV interfero- 
metric imaging framework developed in this work. We attempt to 
recover the sky intensity directly on the sphere x s e C Ns by solving 
the inverse problem 

j = O s x s + n, (11) 

with 

O s = WMFCAP. 

The measurement operator on the sphere O s e C MxNs simply con- 
sists of augmenting the operator on the plane O p by prepending a 
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projection from the sphere to the plane P, which incorporates a con- 
volutional re-gridding as discussed in Sec. |3.2| Careful considera- 
tion is given to samplings on the sphere and plane to ensure that the 
planar grid is sampled sufficiently to accurately represent the pro- 
jection of a band- limited signal defined on the sphere (as discussed 
in Sec. |3.1| >. For incomplete visibility coverage the inverse problem 
defined by |TTJ becomes ill-posed. We solve this ill-posed inverse 
problem in a compressed sensing framework by solving the BP and 
TV minimisation problems defined by {6} and {7} respectively. In 
this work we consider the Dirac sparsity basis on the sphere for BP 
reconstructions. 




(a) Real part 



(b) Imaginary part 



We expect the performance of compressed sensing reconstruc- 
tions to be enhanced by recovering the sky intensity in the space 
where it naturally lives {i.e. on the sphere), since we expect spar- 
sity to be reduced in this space. Furthermore, the effectiveness of 
the spread spectrum phenomenon is enhanced by going to a WFOV. 
As discussed in Sec. |2.3[ the greater the frequency content of the 
w-term, the more effective the spread spectrum phenomenon. By 
examining the w-term plotted in Fig. [5] we see that its frequency 
content increases with distance from the origin. Consequently, the 
larger the FOV, the higher the frequency content achieved by the 
w-term and the more effective the spread spectrum phenomenon. 
Moreover, a secondary enhancement arises by eliminating the first 
order assumption \\l\\* w «: 1 made previously by Wiaux et al. 
(2009b), since this assumption reduces the maximum frequency 
content of the w-term achieved on a given FOV (the band-limits 
of the various w-terms are defined explicitly in Sec. |4.1| ). Finally, 
we comment on the impact of mutual coherence on reconstruction 
performance. Since we consider the Dirac sparsity basis, coherence 
is optimal on the plane (recall that the Dirac and Fourier bases 
are maximally incoherent). However, the projection of Dirac ba- 
sis functions defined on the sphere to the plane does not result in 
Dirac functions on the plane (due to the convolutional re-gridding). 
Consequently, coherence in the spherical setting is suboptimal. The 
spread spectrum phenomenon acts to increase incoherence only in 
the case where it is not already maximally incoherence, thus we ex- 
pect the spread spectrum phenomenon to be ineffective in the planar 
setting but to improve reconstruction performance in the spherical 
setting. For TV reconstructions, although a sparsity basis does not 
exist, one may gain some intuition regarding the impact of coher- 
ence by the following argument. For a piecewise constant signal 
that is sparse in the magnitude of its gradient (i.e. has a gradient 
defined by Dirac functions), the spectrum of the magnitude of its 
gradient must be flat. The gradient operator in space essentially cor- 
responds to a multiplication by frequency in Fourier space, hence 
the original spectrum of the piecewise constant signal must evolve 
as the inverse of frequency. Since this differs to the optimally inco- 
herent spectrum which is flat, the coherence must be suboptimal for 
TV reconstructions. We therefore expect the spread spectrum phe- 
nomenon to provide improvements both on the sphere and plane 
for TV reconstructions, with a greater enhancement expected on 
the sphere due to the spatial spreading of the projection operator. In 
any case, sparsity typically has a greater impact on the performance 
of compressed sensing reconstructions than coherence. These con- 
siderations pertaining to sparsity, mutual coherence and the spread 
spectrum phenomenon lead us to expect an improvement in the per- 
formance of compressed sensing reconstructions when recovering 
interferometric images in the WFOV framework developed here. 



Figure 3. Real and imagin ary p arts of the w-term modulation C(||/||2) (for 
Wd = 1/ V2; see text Sec. |4.1) . Notice that the frequency content of the 
w-term modulation increases with distance from the origin (image centre). 
Consequently, for a WFOV the w-term modulation spreads the spectrum 
more effectively, enhancing the performance of the spread spectrum phe- 
nomenon. Dark and light regions correspond to positive and negative values 
respectively. 



4 SIMULATED RECONSTRUCTIONS 

We evaluate the performance of the WFOV interferometric imaging 
framework defined on the sphere, as outlined in Sec. [3] making 
a direct comparison with planar reconstructions. After describing 
the observational set-up, performance is quantified thoroughly in 
a low-resolution setting on sets of simulations of sources with a 
Gaussian profile. A more realistic setting at a higher resolution is 
then considered, where reconstruction performance is evaluated on 
a single simulated observation of Galactic dust emission. 



4.1 Observational set-up 

Simulated observations of real signals are made on the FOV de- 
fined by the angular opening # FO v = 90°, corresponding to a planar 
FOV of L = V2. A real Gaussian primary beam is assumed, with 
full-width-half-maximum (FWHM) of one half of the field of view. 
Random visibility coverage is considered, with visibility measure- 
ments falling on the discrete planar grid of spatial frequencies u t . 
We consider incomplete visibility coverage, with only 2-25% of 
the discrete visibilities measured (since we consider real signals, 
measuring 50% of the complex visibilities corresponds to a num- 
ber of measurements identical to the number of unknowns in the 
real signal that we attempt to recover). 

As discussed previously, although we consider a WFOV we 
restrict ourselves to the case of constant w. We parameterise the 
continuous w in terms of the discrete compon ent w d , following the 



Wiaux et al. 



( 2009a| ) ofw = w d N p 1/2 /L 2 . 



parameterisation used by 

The band-limit of the modulating w-term C(||/|| 2 ) is given approxi- 
mately by 



f r(w) 



2LVRW 



corresponding to its maximum instantaneous frequency. Under the 
first order small FOV assumption Hfll* w« 1, the w-term reduces 
to the linear chirp Ci(||/|| 2 ) with approximate band-limit 

cf> w d N p l/2 



2L 

Notice that the band-limit for given non-zero values of w d , N p and 
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L is greater in the absence of the small FOV assumption, justify- 
ing rigorously the secondary enhancement discussed in Sec. |3.4| 
of the spread spectrum phenomenon due to the WFOV. The band- 
limit of the spread signal is given by the sum of the original band- 
limit w max = N v l/2 /2L and the band-limit of the w-term. Previously 
[Wiaux et aL] ( |2009b"t considered the linear chirp Ci(p| 2 ) and the 
value Wd = 1, corresponding to spreading by a factor of two. We 
also consider spreading by a factor of two, but since we consider 
the exact w-term C(||/|| 2 ), this corresponds to a value of w d = 1/ V2. 
The corresponding continuous w is of the same order as the maxi- 
mum visibility measurements in u and v, i.e. w = w max , hence it is 
an appropriate value to consider when studying the spread spectrum 
phenomenon. Note that in the absence of a WFOV, the value of w 
considered by Wiaux et al. (2009b) to achieve the same spreading 
was a factor of 2/L greater than w max . Since the band-limit of the 
original signal is doubled due to modulation by the w-term, to avoid 
aliasing we apply an upsampling operator to increase the resolution 
of the planar grid by a factor of two prior to application of the mod- 
ulation (upsampling is performing by zero-padding in the Fourier 
domain). In the subsequent analysis we consider the exact w-term 
Cfllflb), with values of w d = {0, 1/ V2} to highlight the effective- 
ness of the spread spectrum phenomenon. 

Instrumental noise is also added to the simulated visibilities. 
Independent identically distributed Gaussian noise is assumed with 
variance cr 2 n = 10~ 3 a 2 , where cr 2 y is the variance of the visibilities in 
the absence of noise and w-term modulation. The added instrumen- 
tal noise results in observed visibilities at a signal-to-noise ratio of 
SNR n = 101og 10 (a^/0 = 30dB. 

4.2 Gaussian sources 

The WFOV interferometric imaging framework is evaluated thor- 
oughly in this section on simulated observations of Gaussian 
sources. These simulations are first described, before we analyse 
their sparsity properties. Reconstruction performance is then eval- 
uated both in the spherical and planar settings. Finally, reconstruc- 
tion performance is compared to calculations of cumulative coher- 
ence. 



4.2.1 Simulations 

In order to perform a thorough evaluation, we analyse sets of simu- 
lations of Gaussian sources of various size, with each set containing 
30 simulations, for all variations of reconstruction procedures (BP 
and TV reconstructions on the plane and the sphere, both with and 
without application of the spread spectrum phenomenon) and for 
various visibility coverages. Due to the large number of simulated 
reconstructions, we restrict these simulations to a low resolution. 
We consider a HEALPix resolution parameter of Af si d e = 32, cor- 
responding to a harmonic band-limit of ^ max = 88. The remaining 
parameters defining the resolution of the simulations then follow 
from the considerations discussed in Sec. |3.1| we find N s = 1740, 
Umax = 19.6 and N p = 58 x 58 = 3364. The size of the Gaussian 
kernel in the convolutional re-gridding of the projection operator, 
defined by its standard deviation o- P , is chosen to ensure that the 
kernel is well sampled. Since a kernel with small support is re- 
quired, a local tangent plane approximation is made to relate its 
Fourier size to its spatial size through <x F = (2^cr P ) _1 . We ensure 
2cr F is within the band-limit supported by the planar grid, resulting 
in the kernel size cr P = 0.02 radians. 

For each simulation, 10 Gaussian sources of the same size 



are laid down at random positions within the FOV and with 
random amplitudes A s e [0.5,1.0]. For each set of 30 simu- 
lations, source objects of a different size are considered, with 
sizes defined by the standard deviation of the Gaussian source 
os e {0.01,0.02, 0.04, 0.10} radians. For some cases the standard 
deviation of the source is smaller than the pixel size, hence the 
Gaussian sources are not necessarily well sampled on the spheri- 
cal grid. However, this is of no concern since the simulated sig- 
nals are defined by their discrete version and no contact is subse- 
quently made with the continuous representation from which they 
originate. To simulate compact Gaussian objects on the sphere we 
use the COMB ( [McEwen et al.|2008) and S2 ( |McEwen et al.|2007) 
package^] Random visibility coverages are considered with only 
k e {2, 4, 6, 10, 15, 20, 25}% of the discrete visibilities measured. 

4.2.2 Sparsity 

We perform tests to determine whether our hypothesis holds that 
sparsity is reduced by going to the space in which the signal in- 
herently lives. Although sparsity levels will depend highly on the 
signal considered, since each set of Gaussian simulations has simi- 
lar properties we expect sparsity levels to be reasonably consistent 
over the set of simulations. We therefore average the spar sides mea- 
sured over the set of 30 simulations for each source size and also 
provide an indication of their spread. 

Since the signals analysed are compressible rather than ex- 
actly sparse, we require a measure of compressibility. We construct 
such a measure as follows. For the BP problem, we first order the 
coefficients of the signal in the sparsity basis, which in the case of 
the Dirac sparsity basis that we adopt corresponds to simply order- 
ing the sampled signal values. We then set the smallest coefficient 
to zero and measure both the sparsity of the resultant signal and 
also the SNR of the original signal relative to the error between the 
original and resultant signal. We set the next smallest coefficient to 
zero and repeat these two measurements, repeating the procedure 
until all signal values are set to zero. Following this approach we 
build a curve of sparsity against SNR. For the TV problem, we re- 
peat the same procedure but in the space of the magnitude of the 
gradient of the signal, rather than a sparsity basis. 

Sparsity measurements are computed for the signal on the 
sphere and its projected version on the plane in order to make com- 
parisons. Curves are plotted in Fig. [4] for the two extreme source 
sizes (curves for intermediate source sizes are essentially interpo- 
lations between these extremes). Sparsity is clearly enhanced on 
the sphere. Moreover, the ratio of sparsities between the plane and 
sphere (K v /K s ) does indeed approximately follow the ratio of the 
number of samples between the plane and sphere (N p /N s ~ 1.9 for 
#fov = 90°), as predicted naively in Sec. |3.1| Although sparsity will 
always be highly dependent on the signal under investigation, we 
have at least demonstrated for Gaussian sources that the sparsity of 
the signal is indeed enhanced on the sphere. 

4.2.3 Reconstruction performance 

We evaluate reconstruction performance in the WFOV setting by 
recovering interferometric images directly on both the sphere and 
the plane, in order to made a direct comparison. Furthermore, we 
consider BP and TV reconstructions, both with and without ap- 
plication of the spread spectrum phenomenon. It is of particular 
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Figure 4. Imposed sparsity of simulations of Gaussian sources as a func- 
tion of the SNR of the original compressible signal compared to the im- 
posed sparse signal (as explained in the text). Curves correspond to the 
mean of 30 simulations, while shaded regions show one standard deviation 
confidence intervals. Curves are plotted for the signal defined on the sphere 
(red/dot-dashed line) and for its projected version on the plane (blue/solid 
line). Sparsity is clearly enhanced on the sphere. 



interest to see whether our predictions are demonstrated through 
reconstruction quality. 

The measurement contraint e is defined by the level p = 0.99 
(as described in Sec. |2.2| ) and we solve the BP and TV minimisation 
problems using the algorithms derived by |Combettes & Pes quet 
( 2007) and Durand et al. (2010) respectively. For the low-resolution 
simulations performed here, the computation time required to solve 
the optimisation problems are typically of the order of one minute 
on a standard laptop with a 2.66GHz Intel Core 2 Duo processor 
with 4GB of memory. 

The quality of reconstruction is measured by the SNR 
of the original signal relative to the difference between 
the original and reconstructed signal, defined explicitly by 
SNR S = 101og 10 (a^ s /cr^ where o~l s is the variance of the 

original signal x s and ^ s _ Xs * is the variance of the discrepancy 
signal x s - x s *. Since the sky intensity to be recovered lives inher- 
ently on the sphere, we consider the SNR defined on the sphere. It 
is therefore necessary to lift the image reconstructed on the plane 
to the sphere in order to make a comparison. We do this through 
a projection operator that is the direct analogue of the operator P 
that projects the sphere to the plane, using an identical kernel in 
the convolutional re-gridding. We then compare SNR S , measured 
on the sphere, for both the spherical and planar based reconstruc- 
tions. Fig. [5] shows reconstruction quality measured by SNR S for 
various visibility coverages, averaged over the 30 simulations for 
each source size. Reconstruction quality in the spherical setting is 
clearly superior to the quality of planar reconstructions. However, 
for sources of small size, lifting the planar reconstruction to the 
sphere introduces error and limits the effectiveness of the recon- 



struction. Before discussing reconstruction performance in more 
detail, we consider the SNR defined on the plane. 

We also examine the SNR defined on the plane by 
SNR D = 10 log 10 (<rl Ig 2 .), where crl and a 2 *arethevari- 
ances of the original and discrepancy signals on the plane respec- 
tively. To compute SNR p for interferometric images recovered on 
the sphere, the spherical reconstructions are projected onto the 
plane by the projection operator P. Fig. [6] shows reconstruction 
quality measured by SNR p for various visibility coverages, aver- 
aged over the 30 simulations for each source size. The superiority 
of reconstructions on the sphere is again clear. Even if one were in- 
terested in planar reconstructions, superior reconstruction quality is 
achieved by first recovering interferometric images on the sphere, 
before projecting the recovered spherical image to the plane. In any 
case, we advocate the direct use of spherical reconstructions since 
signal content is not distorted by any projection. 

The reconstruction performance observed in Fig.|5]and Fig. [6] 
is now discussed in more detail and related to the predictions that 
we made through intuitive reasoning. For BP reconstructions, the 
improvement in spherical reconstruction quality due to the spread 
spectrum phenomenon is apparent, whereas the spread spectrum 
phenomenon is clearly ineffective on the plane. For TV reconstruc- 
tions, the spread spectrum phenomenon is effective both on the 
sphere and the plane, although the enhancement on the sphere is 
slightly larger than that on the plane. Also note that the variance of 
reconstruction quality (as indicated by the size of the error bars) is 
reduced by the spread spectrum phenomenon in the cases where it 
is effective. This effect was also observed by Wiaux et al. ( 2009b). 
Notice that the performance of BP reconstructions drops as the size 
of the Gaussian sources increase due to the corresponding increase 
in sparsity value. However, this is not the case for TV reconstruc- 
tions. Although TV reconstructions do not perform as well for sig- 
nals that are extremely sparse in the Dirac basis, they are more sta- 
ble as sparsity varies and are certainly superior for reconstructing 
diffuse signal content. In summary, all of our intuitive predictions 
are indeed manifest in the reconstruction quality observed and the 
superiority of reconstructing the sky intensity directly on the sphere 
in the WFOV interferometric imaging framework that we develop 
is clear. 



4.2.4 Comparison of performance with sparsity and coherence 

Although mutual coherence is a useful measure for theoretical and 
intuitive considerations, it is not well suited to numerical compu- 
tation (as discussed in Sec. |2.2| ). Cumulative coherence is a more 
robust measure of coherence and thus is more suitable for numeri- 
cal evaluation. Furthermore, cumulative coherence is signal depen- 
dent and incorporates sparsity information. We therefore use cumu- 
lative coherence to study the combined sparsity and coherence of 
the interferometric imaging framework in the context of Gaussian 
sources and relate this to the reconstruction performance presented 
in the previous section. However, we may only study cumulative 
coherence in the presence of a sparsity basis, thus the analysis and 
discussion here is restricted to BP reconstructions only. 

In a strict compressed sensing framework with orthogonal 
measurement and sparsity bases and in the absence of noise, the 
number of measurements required to reconstruct a signal accurately 
from incomplete random measurements evolves as {8}. The square 
of the normalised cumulative coherence v ^^N provides an approx- 
imate relative measure of the number of measurements required 
to recover a signal, or similarly, the quality of the reconstruction 



10 J.D. McEwen & Y. Wiaux 





(j) Typical simulation (<xs = 0.10) (k) BP reconstruction (<xs = 0.10) (1) TV reconstruction (<xs = 0.10) 



Figure 5. Reconstruction quality for simulated Gaussian sources measured by the SNR defined on the sphere (SNR S ). Each row of panels corresponds to 
Gaussian sources of a given size. The first column of panels shows a typical simulation of Gaussian sources on the sphere. The remaining columns of panels 
show reconstruction quality, with the second column illustrating the performance of BP reconstructions, while the third column illustrates the performance of 
TV reconstructions. Curves are plotted for reconstructions performed on the sphere (red/diamonds) and on the plane (blue/circles), both in the absence (solid 
lines) and presence (dashed lines) of the spread spectrum phenomenon. Reconstruction quality is averaged over 30 simulations for each source size and error 
bars corresponding to one standard deviation are shown. 
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(a) Typical simulation (cr s = 0.01) (b) BP reconstruction (<xs = 0.01) (c) TV reconstruction (cr s = 0.01) 




(d) Typical simulation (cr s = 0.02) (e) BP reconstruction (<xs = 0.02) (f) TV reconstruction (cr s = 0.02) 
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(g) Typical simulation (<ts = 0.04) 



(h) BP reconstruction (cr s = 0.04) 



(i) TV reconstruction (crs = 0.04) 



(j) Typical simulation (crs =0.10) 
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(k) BP reconstruction (crs =0.10) 
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(1) TV reconstruction (crs =0.10) 



Figure 6. Reconstruction quality for simulated Gaussian sources measured by the SNR defined on the plane (SNR p ). Each row of panels corresponds to 
Gaussian sources of a given size. The first column of panels shows a typical simulation of Gaussian sources projected onto the plane. The remaining columns 
of panels show reconstruction quality, with the second column illustrating the performance of BP reconstructions, while the third column illustrates the 
performance of TV reconstructions. Curves are plotted for reconstructions performed on the sphere (red/diamonds) and on the plane (blue/circles), both in the 
absence (solid lines) and presence (dashed lines) of the spread spectrum phenomenon. Reconstruction quality is averaged over 30 simulations for each source 
size and error bars corresponding to one standard deviation are shown. 
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performance for a given number of measurements. Although the 
interferometric imaging framework we consider differs from the 
theoretical compressed sensing framework, the normalised cumu- 
lative coherence nevertheless provides a measure of the impact of 
sparsity and coherence on reconstruction performance. 

We compute cumulative coherences for all of the sets of Gaus- 
sian simulations, averaging over the 30 simulations for each source 
size. Normalised cumulative coherences are computed on the plane 
and sphere, denoted v p (w d ) yjN^ and v s (w d ) V^s respectively, both 
in the presence (w d = 1 / V2) and absence (w d = 0) of the spread 
spectrum phenomenon. Results are displayed in Table [T] Four in- 
sights can be made by studying the coherence values reported in 
Table[T] Firstly, the coherences on the plane are not significantly al- 
tered in the presence of the spread spectrum phenomenon, indicat- 
ing that the spread spectrum phenomenon is likely to be ineffective 
in improving reconstruction performance in this setting. Secondly, 
the coherences on the sphere are reduced by application of the 
spread spectrum phenomenon, indicating that the spread spectrum 
phenomenon is likely to be effective in this setting. Thirdly, we note 
that coherences on the sphere are always lower than the correspond- 
ing values on the plane, highlighting the superiority of spherical re- 
construction. Fourthly, we note that cumulative coherences increase 
with source size, indicating that reconstruction performance should 
reduce. All of these findings verify our intuitive expectations and 
are consistent with the reconstruction performance presented in the 
previous section. 

4.3 Galactic dust 

Now that the WFOV interferometric imaging framework developed 
in this article has been evaluated thoroughly, in this section we sim- 
ply illustrate the reconstruction, at a higher resolution, of a more re- 
alistic simulation of Galactic dust emission. The simulation is first 
described, before reconstruction performance is evaluated. 

4.3.1 Simulation 

For this higher resolution simulation we consider the 94GHz map 
of predicted submillimeter and microwave emission of diffuse in- 
terstellar Galactic dust (Finkbeiner et al. 1999 ), hereafter referred 
to as the FDS map. This predicted map is based on the merged 
Infrared Astronomy Satellite (IRAS) and Cosmic Background Ex- 
plorer Diffuse Infrared Background Experiment (COBE-DIRBE) 
observations produced by Schlege l et al.| {1998 ). An undersam- 
pled version of the FDS map is available from the Legacy Archive 
for Microwave Background Data Analysis (LAMBDA) in the 
HEALPix pixelisation at resolution A^ side = 512. 

We consider the same observational set-up discussed in 
Sec. |4.1[ focusing on the FOV centred on Galactic coordinates 
(l,b) = (210°, -20°). The full-sky FDS map and the FOV con- 
sidered are illustrated in Fig. [7] We downgrade the original FDS 
map to a resolution of A^ side = 128 for our simulated observations, 
corresponding to a harmonic band-limit of £ max = 348. All other 
resolution parameters follow once the harmonic band-limit on the 
sphere and the size of the FOV are chosen; we find N $ = 28560, 
"max = 78.4 and N p = 214 x 214 = 45796. The size of the Gaussian 
kernel in the convolutional re-gridding of the projection operation 
is chosen by the same condition as before, resulting in the value 



9 http://lambda.gsfc.nasa.gov/ 




(a) Mollweide projection of full- sky 




(b) Orthographic projection of FOV 



Figure 7. FDS map of predicted submillimeter and microwave emission of 
diffuse interstellar Galactic dust. The FOV of angular opening #fov = 90° 
and centred on Galactic coordinates (/, b) = (210°, -20°) is considered for 
simulating interferometric observations. In panel (a) the FOV is located near 
the Galactic plane towards the left edge of the image. 

cr P = 0.004 radians. Random visibility coverage is again consid- 
ered, with only 25% of the discrete visibilities measured. Note that 
the incomplete visibility coverage creates a very challenging setting 
for the recovery of a realistic, diffuse image. 

4.3.2 Reconstruction performance 

All of the reconstruction techniques discussed previously are ap- 
plied to reconstruct the FDS map on the sphere (i.e. BP and TV 
reconstructions in the spherical and planar settings, both with and 
without application of the spread spectrum phenomenon). The 
same optimisation algorithms and measurement constraint level p 
considered previously are applied. For these high-resolution sim- 
ulations, the computation time required to solve the optimisation 
problems are typically of the order of 15 minutes on the same ma- 
chine described previously. The results presented in Sec. [42] indi- 
cate that for diffuse signals the SNR measured on the sphere (SNR S ) 
is not adversely affected significantly by lifting planar reconstruc- 
tions to the sphere. Due to the diffuse nature of the underlying sig- 
nal (and also since it is inherently defined on the sphere), we there- 
fore measure reconstruction quality by SNR S only^Note that for 
such signals the SNR metric is not a perfect measure of reconstruc- 
tion quality and a visual inspection remains important. The diffuse 
nature of the FDS map and the results presented previously suggest 

10 Although we only quote SNR measured on the sphere, the SNR mea- 
sured on the plane (SNR p ) were also examined. Corresponding SNR values 
do change marginally, nevertheless the conclusions drawn remain the same 
regardless of whether SNR S or SNR p is examined. 



Compressed sensing for wide-field imaging 13 



Source 




Plane 






Sphere 




Difference 












v s (w d ) VA^ 




A(w d ) 




w'd = 


Wd= Tl 


A P 


w'd = 


Wd = Ti 


A s 


w d = w d = ±= 


0.01 


8.14 + 0.72 


8.31 ±0.78 


0.17 ±0.10 


5.07 ± 0.52 


4.69 ± 0.55 


-0.38 ± 0.06 


-3.08 ± 0.23 -3.62 ± 0.25 


0.02 


11.05 ±0.95 
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12.84 ± 1.25 


-1.01 ±0.11 


-2.46 ±0.07 -3.52 ±0.13 


0.10 


22.32 ± 0.25 


22.33 ± 0.25 


0.01 ±0.01 


19.79 ± 0.25 


18.23 ±0.25 


-1.73 ±0.03 


-2.53 ±0.02 -4.12 ±0.05 



Table 1. Normalised cumulative coherences, on the plane and sphere, of simulations of Gaussian sources for various size os- Coherences values computed 
both in the presence (wd = 1 / V2) and absence (wd = 0) of the spread spectrum phenomenon are displayed. Normalised coherence differences between the 
presence and absence of the spread spectrum phenomenon are computed on the plane and sphere (with A p and A s denoting, respectively, the coherence for 
Wd = subtracted from the coherence for w d = 1/ V2). In addition, normalised coherence differences between the plane and sphere are also computed, with 
A(wd) = Vp(wd) JN^ - v s (wd) ^^N^. Quoted values correspond to the mean of 30 simulations, with errors corresponding to one standard deviation. Notice that 
coherence differences are more stable than coherence values over the simulations. These measures incorporate both the impact of sparsity and coherence on 
reconstruction performance (with lower values indicating superior performance). 



that the TV reconstruction on the sphere in the presence of spread 
spectrum phenomenon will be most effective. 

Reconstructed spherical interferometric images of the FDS 
map are displayed in Fig. [8] with SNR S of each recovery also spec- 
ified. The quality of BP reconstructions is relatively poor (as ex- 
pected since the FDS map is not particularly sparse in the Dirac 
basis), and a significant difference between planar and spherical 
reconstructions is not readily apparent. Spherical reconstructions 
show greater detail than the planar reconstructions, particularly 
within the centre of the FOV, but suffer near the extremes of the 
FOV since the signal is poorly constrained here due to the low mag- 
nitude of the primary beam. Planar reconstructions also suffer from 
this effect, however it is not as apparent in reconstructed spherical 
images since a small degree of smoothing is performed in the pro- 
jection operator when lifting the reconstructed planar signal to the 
sphere. The quality of reconstruction is improved considerably for 
the TV case, again as expected since the diffuse FDS map is much 
sparser in the magnitude of its gradient than it is in the Dirac basis. 
Spherical reconstructions again show greater detail than the pla- 
nar reconstructions, however in the absence of the spread spectrum 
phenomenon SNR S is in fact lower for the spherical reconstruc- 
tion. On visual inspection the spherical reconstruction is clearly 
superior as structure of much finer detail is reconstructed, high- 
lighting the weakness of the SNR metric. In any case, the superi- 
ority of the spherical reconstruction is clearly apparent for the TV 
reconstruction when the spread spectrum phenomenon is present, 
both through visual inspection and comparison of SNR S . This is 
the most effective reconstruction technique for the FDS map for 
both planar and spherical reconstructions, as expected. Comparing 
the most effective reconstruction on the plane and sphere, recover- 
ing the sky intensity directly on the sphere improves reconstruction 
quality from 13.7dB to 19.3dB. 



5 CONCLUSIONS 

Incorporating WFOV contributions when reconstructing interfero- 
metric images is becoming increasingly important, particularly for 
next-generation interferometers, such as the SKA and upcoming 
pathfinder telescopes. If these contributions are ignored, the fidelity 
of reconstructed images may suffer. In this article we have extended 
the compressed sensing interferometric imaging techniques devel- 
oped by |Wiaux et aT] ( |2009a| ) and |Wiaux et"aT]p009b| > to a WFOV 
In doing so, we recover interferometric images defined directly on 



the celestial sphere. We augment the usual measurement operator 
with a projection from the sphere to the plane, which essentially 
corresponds to a change from Cartesian to spherical coordinates. 
Practically, however, the projection is complicated by the discrete 
setting and careful consideration is given to the sampling of signals 
on the sphere and plane to ensure that the planar grid is sufficiently 
sampled to support the maximum projected frequency content of a 
band-limited signal on the sphere. Although a projection is incor- 
porated in this framework, it is included in the measurement opera- 
tor and thus is regularised when solving the interferometric inverse 
problem. Moreover, it is the space where recovery is performed 
that drives the performance of compressed sensing reconstructions, 
through sparsity and coherence considerations. The framework re- 
mains general and does not rely on any small FOV assumptions. 

The effectiveness of the spread spectrum phenomenon is en- 
hanced when going to a WFOV, while sparsity is promoted by 
recovering images directly on the sphere. These predictions have 
been verified by numerical tests and are also manifest in recon- 
struction performance. Low-resolution simulations of Gaussian 
sources were considered to quantify reconstruction performance 
thoroughly. Interferometric images were recovered directly on the 
sphere and the plane in order to make comparisons. For simulated 
images extremely sparse in the Dirac basis, BP reconstructions 
were shown to perform very well. However, as Dirac sparsity was 
reduced the quality of BP reconstructions fall, while the quality of 
TV reconstructions remained relatively stable. For diffuse images, 
TV reconstructions were shown to be superior since such signals 
are much more sparse in the magnitude of their gradient than in 
the Dirac basis. In all cases, the superior quality of recovering in- 
terferometric images directly on the sphere was clear. A simulation 
of diffuse interstellar Galactic dust was then performed to demon- 
strate WFOV reconstruction techniques in a more realistic, higher 
resolution setting. For the diffuse FDS map considered, TV recon- 
struction on the sphere in the presence of the spread spectrum phe- 
nomenon was most effective, as expected. For this case, reconstruc- 
tion quality was improved from 13.7dB for the planar reconstruc- 
tion to 19.3dB when recovering the interferometric image directly 
on the sphere. 

The compressed sensing techniques developed for interfero- 
metric imaging by Wiaux et al. ( 2009a) and Wiaux et al. ( 2009b), 
and extended here to a WFOV, remain somewhat idealised. Ran- 
dom visibility coverage is assumed, with the spatial frequencies 
probed by the interferometer also assumed to fall on discrete grid 
points. Furthermore, to study the spread spectrum phenomenon a 
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Figure 8. Reconstructed spherical interferometric images of the FDS map. The first column of panels shows images of planar recoveries lifted to the sphere, 
while the second column shows images recovered on the sphere directly. Reconstructions both in the presence (w d = 1 / V2) and absence (w d = 0) of the 
spread spectrum phenomenon are displayed. 
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constant w is assumed. These restrictions have been necessary to 
remain as close to the theory of compressed sensing as possible 
during the development and evaluation of interferometric imaging 
techniques. In reality, however, w will not be constant and the per- 
formance of the spread spectrum phenomenon will lie between the 
extreme cases that we have considered of w = and w = u max . 
Extensions to realistic and continuous visibility coverage and their 
impact on compressed sensing based interferometric imaging are 
now of considerable importance. In general, compressed sensing 
addresses imaging by optimising both reconstruction and acquisi- 
tion, while we have essentially focused on reconstruction only. The 
possibility of optimising the configuration of interferometers to en- 
hance the spread spectrum phenomenon for compressed sensing 
reconstruction is an exciting avenue of research at the level of ac- 
quisition. In addition, direction dependent beam effects may also 
provide an alternative source of the spread spectrum phenomenon. 
All of these issues should be studied in future works. Furthermore, 
the performance of other sparsity bases on the sphere, such as 
scale discretised wavelets (Wi aux et al.|20 08 ), should also be stud- 
ied. Next-generation radio interferometers will inherently observe 
very large fields of view, thus WFOV interferometric imaging tech- 
niques, such as the compressed sensing techniques developed in 
this article and the future research outlined here, are of increasing 
importance to ensure that the fidelity of reconstructed images keeps 
pace with the capabilities of new instruments. 
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